#' Univariate and Bivariate Mark correlation functions and random labeling test
#' 
#' @param focus the community object contains focus individuals
#' @param neighbors the community object contains all neighors, if NULL the
#'   function become univariate mark correlation funciton, otherwise bivariate
#'   case.
#' @param r the spatial distance to calcualte the mark function
#' @param h bandwidth used for smoothing the mark correlation function
#' @param testfun name of different test functions, see details below
#' @param isnormal a logical flag to indicate whether the mark correlation
#'   funciton should be normalized
#' @param edgeCor the method of edge correction, currently only 'border'
#'   correction was implemented
#' @param nullmodel_fun the name of null model function, currently 'random_mark'
#'   null model was implemented, see details below
#' @param random_focus if it is true, the null model will apply on focus,
#'   otherwise on neighbors
#' @param nrep number of replicates in the null model
#' @param alpha,alternative the significant level and ways in calculate the
#'   confidence interval, see details in conf
#' @param savenull a logical flag to indicate whether the values of null model
#'   should be returned.
#' @param mc.cores number of cores used. default is 1. it only works under linux OS.
#' @param ... only for advanced user, other parameters that can be passed to
#'   use defined testfun, meanfun and nullmodelfun
#'   
#' @details There is an general interesting to find out whether the marks of
#' points show some sort of spatial correlation, conditional on the locations of
#' the corresponding unmaked pattern. This is a full package of calculating mark
#' correlation function with various test fuctions and various null models.
#' 
#' Currently, six test functions have been implemented. They are 
#' \tabular{ll}{ 
#' 1. \tab \eqn{t(m_i,m_j)=m_i*m_j} \cr 
#' 2. \tab \eqn{ t(m_i,m_j)=mi }, r-mark function \cr 
#' 3. \tab \eqn{t(m_i,m_j)=mj} \cr 
#' 4. \tab \eqn{t(m_i,m_j)=(m_i-m_j)^2/2}, mark variogram, please normalize marks for each type for meaningful comparison \cr 
#' 5. \tab \eqn{t(m_i,m_j)=(m_i-mu)(m_j-mu)} \cr 
#' 6. \tab \eqn{t(m_i,m_j,r)=(m_i-mu(r))(m_j-mu(r))} \cr
#' 7. \tab \eqn{t(m_i,lambda,g_i(r))=(m_i-mu)(lambda*g_i(r)-lambda*g(r)) },  only consider the density of neighbors, nonaccumulative version \cr 
#' 8. \tab \eqn{t(m_i,lambda,k_i(r),dbh_j)=(m_i-mu)(lmabda*k_i(r)*dbhj-somemean)}, add numeric mark info with neighbor density \cr 
#' 9. \tab \eqn{t(m_i,lambda,k_i(r))=(m_i-mu)(lambda*k_i(r)-lambda*k(r))},  only consider the density of neighbors, accumulative version \cr 
#' }
#' 
#' 
#' @return It will return a 'mcr' object contains four columns: r (spatial
#' scales), and obs, lower and upper values of the mark correlation function it
#' also has an attribute saved all normalized constant in the observed and all
#' simulated patterns in the null model
#' 
#' @author Heming Liu, Guochun Shen
#' 
#' @examples
#' 
#' data(BCI)
#' set.seed(1)
#' BCI$dbh=rnorm(total_abundance(BCI),10,1)
#' focus=get_populations(BCI,"ACALDI")
#' r=seq(0,10,0.5)
#' 
#' #univariate mark correlation function
#' re1=markcorr(focus,r=r,testfun="t1")
#' plot(re1)
#' re2=markcorr(focus,r=r,testfun="t2")
#' plot(re2)
#' #you can try other test functions
#' 
#' #bivariate mark correlation function
#' 
#' neighbors=get_populations(BCI,"BEILPE")
#' re1=markcorr(focus,neighbors,r=r,testfun="t1")
#' plot(re1)
#' 
#' re2=markcorr(focus,neighbors,r=r,testfun="t2")
#' plot(re2)
#' 
#' #you can try other test functions
#' 
#' re8=markcorr(focus,neighbors,r=seq(0,10,0.5),h=0.5,testfun="t8",edgeCor="border",nullmodel_fun="random_mark",nrep=19)
#' plot(re8)
#'
#'
#' #global test for 2m to 5m distance range
#' re7=markcorr(focus,neighbors,r=seq(0,10,0.5),h=0.5,testfun="t7",edgeCor="border",nullmodel_fun="random_mark",nrep=199,savenull=T)
#' mt=cbind(re7$obs,attr(re7,"null_value"))
#' sel=re7$r>=2 & re7$r<=5
#' lmGofTest(mt[sel,]) #return a pvalue
#' 
#' 
#bivariate mark correlation function
#user interface
#' @export
markcorr=function(focus,neighbors=NULL,fmarkname="dbh",nmarkname="dbh",r=seq(0,10,1),h=2,testfun="t1",isnormal=TRUE,
                  edgeCor="border",nullmodel_fun="random_mark",random_focus=TRUE,nrep=199,alpha=0.05,alternative="both",
                  savenull=FALSE,mc.cores=1,...){
  
  if(is.null(neighbors)){
    biver=FALSE
    neighbors=focus
    #if(fmarkname!=nmarkname)
    #  stop("the mark name should be the same in univariate case")
  }else{
    biver=TRUE
    
  }
  
  #check type of data 
  if(!(("community" %in% class(focus)) | ("population" %in% class(focus))) ) 
    stop("class of focus/neighbors is not community/population")
  if(!(("community" %in% class(neighbors)) | ("population" %in% class(neighbors))) ) 
    stop("class of focus/neighbors is not community/population")
  
  #check the name of marks
  if(!(fmarkname %in% get_trait_names(focus)) | !(nmarkname %in% get_trait_names(neighbors)))
    stop("can't find the given marks in the focus and neighbors")
  
  #total number of focus individuals
  N=total_abundance(focus)
  #minimum distance to border, prepare for border correction
  min_bdist=dist_to_border(focus)
  #total number of distance bins
  nr=length(r)
  
  #adjust r
  rmax=r[nr]+h
  rmin=r[1]-h
  if(rmin<0) rmin=0
  #calculate spatial neighbors
  fr=frnn(neighbors,focus$x,focus$y,c(rmin,rmax))
  
  #devide neighbors into different distance bins
  def_cut=define_rectSmooth(r,h)
  fr_cut=lapply(fr,function(x) apply_rectSmooth(def_cut,x$dist))
  #change the distance index to the index of individuals in the neighbors object
  fr_cut=lapply(1:N,function(x) lapply(fr_cut[[x]], function(y) fr[[x]]$index[y]))
  
  #chose right test function, mean value and normalized constant
  test_pars=switch (testfun,
                    "t1" = list(tfun=testfun_t1,ct=mark_mean(focus,fmarkname)*mark_mean(neighbors,nmarkname)),
                    "t2" = list(tfun=testfun_t2,ct=mark_mean(focus,fmarkname)),
                    "t3" = list(tfun=testfun_t3,ct=mark_mean(neighbors,fmarkname)),
                    "t4" = list(tfun=testfun_t4,ct=mark_covariance(focus,neighbors,fmarkname,nmarkname)),
                    "t5" = list(tfun=testfun_t5,ct=mark_covariance(focus,neighbors,fmarkname,nmarkname),
                                mur=c(mark_mean(focus,fmarkname),mark_mean(neighbors,nmarkname))),
                    "t6" = list(tfun=testfun_t6,ct=mark_covariance(focus,neighbors,fmarkname,nmarkname),
                                mur=mark_mean_r(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist)),
                    "t7" = list(tfun=testfun_t7,ct=ind_g_std(focus,fmarkname,fr_cut,r,min_bdist,...),
                                mur=meanfun_t7(focus,fmarkname,fr_cut,r,min_bdist)),
                    "t9" = list(tfun=testfun_t7,ct=ind_K_std(focus,fmarkname,fr,r,min_bdist,...),
                                mur=meanfun_t9(focus,fmarkname,fr,r,min_bdist)),
                    "t8" =list(tfun=testfun_t7,ct=ind_Km_std(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist,...),
                               mur=meanfun_t8(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist))
  )
  
  #select proper null models
  nmfun=switch (nullmodel_fun,
                random_mark = random_mark
  )
  
  #observed mark correlation function
  re_obs=mark_core(focus, neighbors, fmarkname,nmarkname,test_pars[["tfun"]],isnormal, test_pars[["mur"]], fr_cut, fr, min_bdist, nr, N, r,...)
  
  if(isnormal){
    re_obs=re_obs/test_pars[["ct"]]
  }
  
  #mark correlation function under null model
  if(nrep!=0){
    require(parallel)
    re_null=mclapply(1:nrep,function(x){
      
      if(biver){
        if(random_focus){
          focus=nmfun(focus,fmarkname,...)
        }else{
          neighbors=nmfun(neighbors,nmarkname,...)
        }
      }else{
        #in univariate case, both focus and neighbors should be randomized
        focus=neighbors=nmfun(focus,fmarkname,...)
      }
      
      #update mur for testfun8
      if((!biver | !random_focus ) & testfun=="t8"){
        test_pars[["mur"]]=meanfun_t8(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist)
      }
      
      re_sim=mark_core(focus, neighbors, fmarkname,nmarkname,test_pars[["tfun"]],isnormal, test_pars[["mur"]], fr_cut, fr, min_bdist, nr, N, r,...)
      
      #update normalized constant for testfun8 if required
      if(isnormal){
        if( (!biver | !random_focus) & testfun=="t8"){
          ct_sim=ind_Km_std(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist,...)
        }else{
          ct_sim=test_pars[["ct"]]
        }
        re_sim=re_sim/ct_sim
      }
      
    },mc.cores=mc.cores)
    
    re_null=unlist(re_null)
    dim(re_null)=c(nr,nrep)
    re_conf=apply(re_null,1,conf,alpha=alpha,alternative=alternative)
    
    re=data.frame(r=r,obs=re_obs,lower=re_conf[1,],upper=re_conf[2,])
  }else{
    savenull=FALSE
    re=data.frame(r=r,obs=re_obs,lower=NA,upper=NA)
  }
  
  attr(re,"normal_constant")=test_pars[["ct"]]
  if(savenull)
    attr(re,"null_value")=re_null
  class(re)=c("mcf",class(re))
  return(re)
}


mark_core=function(focus, neighbors, fmarkname,nmarkname,testfun,isnormal, mur, fr_cut, fr, min_bdist, nr, N, r,...) {
  
  re=numeric()
  #for each distance bin
  for(ri in 1:nr){
    #initialize numeritor and denumertor for each distance bin
    cri=dei=0
    #for each focus individual
    for(i in 1:N){
      #border correction
      if(min_bdist[i]<r[ri] | length(fr_cut[[i]][[ri]])==0) next()
      #print(testfun(focus,neighbors,fmarkname,nmarkname,i,fr_cut[[i]][[ri]],mur,ri,...))
      cri=cri+testfun(focus,neighbors,fmarkname,nmarkname,i,fr_cut[[i]][[ri]],mur,ri,...)
      dei=dei+length(fr_cut[[i]][[ri]])
    }
    #if(ri==25) browser()
    re[ri]=cri/dei
    
  }
  #re[is.na(re)]=0
  return(re)
}

#t(mi,mj,r)=mi*mj
testfun_t1=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum(get_traits(focus,fmarkname,drop=TRUE)[findex]*get_traits(neighbors,nmarkname,drop=TRUE)[nindex])
}


testfun_t2=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  get_traits(focus,fmarkname,drop=TRUE)[findex]*length(nindex)
}

#t(mi,mj,r)=mj
testfun_t3=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum(get_traits(neighbors,nmarkname,drop=TRUE)[nindex])
}

testfun_t4=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum((get_traits(focus,fmarkname,drop=TRUE)[findex]-get_traits(neighbors,nmarkname,drop=TRUE)[nindex])^2)/2
}

################
testfun_t5=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum((get_traits(focus,fmarkname,drop=TRUE)[findex]-mur[1])*(get_traits(neighbors,nmarkname,drop=TRUE)[nindex]-mur[2]))
}

#####################################
#t(mi,mj,r)=(mi-mu_i(r))(mj-mu_j(r))
testfun_t6=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum((get_traits(focus,fmarkname,drop=TRUE)[findex]-mur[ri,1])*(get_traits(neighbors,nmarkname,drop=TRUE)[nindex]-mur[ri,2]))
}

#####################
#t=(m_i-mu)(lambda*K_i(r)-lambda*K(r))
testfun_t7=function(focus,neighbors,fmarkname,nmarkname,findex,nindex,mur,ri,...){
  sum((get_traits(focus,fmarkname,drop=TRUE)[findex]-mur[[1]])*(mur[[2]][ri,findex]-mur[[3]][ri]))
}


#########################

mark_mean=function(data,markname){
  mean(get_traits(data,markname,drop = TRUE))
}

mark_covariance=function(focus,neighbors,fmarkname,nmarkname){
  sd(get_traits(focus,fmarkname,drop=TRUE))*sd(get_traits(neighbors,nmarkname,drop=TRUE))
}

#mu_i(r)
mark_mean_r=function(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist,...){
  re1=re=numeric()
  N=total_abundance(focus)
  ftrait=get_traits(focus,fmarkname,drop=TRUE)
  ntrait=get_traits(neighbors,nmarkname,drop=TRUE)
  for(ri in 1:length(r)){
    cri=dei=0
    cri1=dei1=0
    for(i in 1:N){
      #border correction
      if(min_bdist[i]<r[ri]) next()
      
      if(length(fr_cut[[i]][[ri]])!=0){
        cri1=cri1+ftrait[i]
        dei1=dei1+1
        tp=ntrait[fr_cut[[i]][[ri]]]
        cri=cri+sum(tp)
        dei=dei+length(tp)
      }
      
    }
    re1[ri]=cri1/dei1
    re[ri]=cri/dei
  }
  return(data.frame(ui=re1,uj=re))
}


meanfun_t7=function(focus,fmarkname,fr_cut,r,min_bdist,...){
  #TODO Consider edge correction
  mu=mean(get_traits(focus,fmarkname,drop=TRUE))
  #non accumulate version
  lk=unlist(lapply(fr_cut,function(x) lapply(1:length(r),function(y) length(unique(unlist(x[y]))) ) ))
  dim(lk)=c(length(r),length(fr_cut))
  for(i in 1:length(min_bdist)){
    lk[r>min_bdist[i],i]=NA
  }
  lkm=apply(lk,1,mean,na.rm=T)
  return(list(mu,lk,lkm))
}

ind_g_std=function(focus,fmarkname,fr_cut,r,min_bdist,...){
  sdi=sd(get_traits(focus,fmarkname,drop=TRUE))
  lk=unlist(lapply(fr_cut,function(x) lapply(1:length(r),function(y) length(unique(unlist(x[y]))) ) ))
  dim(lk)=c(length(r),length(fr_cut))
  for(i in 1:length(min_bdist)){
    lk[r>min_bdist[i],i]=NA
  }
  #sdj=apply(lk,1,sd,na.rm=T)
  sdj=sd(as.numeric(lk),na.rm=T)
  return(sdi*sdj)
}

##
meanfun_t8=function(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist,...){
  #TODO Consider edge correction
  mu=mean(c(get_traits(focus,fmarkname,drop=TRUE)))
  
  mk=get_traits(neighbors,nmarkname,drop=TRUE)
  
  lk=unlist(lapply(fr_cut,function(x) cumsum(unlist(lapply(x,function(y) sum(mk[y]^2/4*pi ))))))
  lk[is.na(lk)]=0
  dim(lk)=c(length(r),length(fr_cut))
  lkm=apply(lk,1,mean,na.rm=T)
  return(list(mu,lk,lkm))
}
#this value will change with null model suffle neighbor traits
ind_Km_std=function(focus,neighbors,fmarkname,nmarkname,fr_cut,r,min_bdist,...){
  sdi=sd(c(get_traits(focus,fmarkname,drop=TRUE)))
  mk=get_traits(neighbors,nmarkname,drop=TRUE)
  
  lk=unlist(lapply(fr_cut,function(x) cumsum(unlist(lapply(x,function(y) sum(mk[y]^2/4*pi ))))))
  sdj=sd(lk,na.rm=T)
  return(sdi*sdj)
}


##
meanfun_t9=function(focus,fmarkname,fr,r,min_bdist,...){
  #TODO Consider edge correction
  mu=mean(c(get_traits(focus,fmarkname,drop=TRUE)))
  
  lk=matrix(nrow=length(r),ncol=length(min_bdist))
  for(i in 1:length(min_bdist)){
    for(ri in 1:length(r)){
      if(min_bdist[i]<r[ri]) next()
      lk[ri,i]=sum(fr[[i]]$dist<ri)
    }
  }
  
  #dim(lk)=c(length(r),length(fr_cut))
  lkm=apply(lk,1,mean,na.rm=T)
  return(list(mu,lk,lkm))
}

ind_K_std=function(focus,fmarkname,fr,r,min_bdist,...){
  sdi=sd(get_traits(focus,fmarkname,drop=TRUE))
  
  lk=matrix(nrow=length(r),ncol=length(min_bdist))
  for(i in 1:length(min_bdist)){
    for(ri in 1:length(r)){
      if(min_bdist[i]<r[ri]) next()
      lk[ri,i]=sum(fr[[i]]$dist<ri)
    }
  }
  
  #dim(lk)=c(length(r),length(fr_cut))
  
  sdj=sd(as.numeric(lk),na.rm=T)
  return(sdi*sdj)
}


#' @export
plot.mcf=function(re,xtext="Spatial distance r",ytext="Mark correlation function"){
  require(ggplot2)
  require(grid)
  conf_data=data.frame(x=c(re$r,rev(re$r)),y=c(re$lower,rev(re$upper)))
  ggplot(data=re)+
    geom_polygon(aes(x=x,y=y),alpha=0.8,fill="grey",data=conf_data)+
    geom_line(aes(x=r,y=obs))+
    ylab(ytext)+xlab(xtext)+
    theme_bw()+
    theme(text=element_text(family = "sans",size = 11),
          plot.margin=unit(rep(0.2,4),"lines"),
          panel.grid.minor=element_blank(),
          panel.grid.major=element_blank())
  
}


markcorr_compare=function(data,rrange=c(1,10),h=0.5,step=NA,edge=c(xmin,xmax,ymin,ymax),index="dbh")
{
  #(rmax-rmin)/2*h should be the integer
  library(scp)
  rmin=rrange[1]
  rmax=rrange[2]
  if(is.na(step)==T)
  {
    int1=(rmax-rmin)/(2*h)
    if(int1%%1!=0)
    {
      print("parameters are not suitable")
    }
  }
  else
  {
    int1=(rmax-rmin)/step
    if(int1%%1!=0)
    {
      print("parameters are not suitable")
    }
  }
  rmin1=rmin+0.000001
  if(rmin-h==0)
  {
    if(is.na(step)==T)
    {
      r1=rmin+2*h
      r2=as.numeric(rmin1)
      while(r1<=rmax)
      {
        r2=c(r2,r1)
        r1=r1+2*h
      }
    }
    else
    {
      r1=rmin+step
      r2=as.numeric(rmin1)
      while(r1<=rmax)
      {
        r2=c(r2,r1)
        r1=r1+step
      }
    }
    
  }
  else
  {
    if(is.na(step)==T)
    {
      r1=rmin+2*h
      r2=as.numeric(rmin)
      while(r1<=rmax)
      {
        r2=c(r2,r1)
        r1=r1+2*h
      }
    }
    else
    {
      r1=rmin+step
      r2=as.numeric(rmin)
      while(r1<=rmax)
      {
        r2=c(r2,r1)
        r1=r1+step
      }
    } 
  }
  data1=data[data$x>=rmax+h+edge[1]&data$x<=edge[2]-rmax-h&data$y>=rmax+h+edge[3]&data$y<=edge[4]-rmax-h,]
  data=community(data$sp,data$x,data$y,plotdim=c(500,400),traits=data.frame(index1=data[[index]]))
  data1=community(data1$sp,data1$x,data1$y,plotdim=c(500,400),traits=data.frame(index1=data1[[index]]))
  n=length(r2)
  cr1=as.numeric()
  for(i in 1:n)
  {
    data.x=data1$x
    data.y=data1$y
    data.m=data1$index1
    data2=frnn(data,data.x,data.y,rRange =c((r2[i]-h),(r2[i]+h)))
    n1=nrow(data1)
    k1=rep(NA,n1)
    k2=k1
    for(j in 1:n1)
    {
      k=data.m[j]*length(data2[[j]]$index)*0.5*h
      k1[j]=k
      k22=length(data2[[j]]$index)*0.5*h
      k2[j]=k22
    }
    if(length(k1[is.na(k1)==T])!=0)
    {
      print("mistake1")
    }
    if(length(k2[is.na(k2)==T])!=0)
    {
      print("mistake2")
    }
    if(sum(as.numeric(k2))==0)
    {
      print("mistake3")
    }
    else
    {
      cr=sum(as.numeric(k1))/sum(as.numeric(k2))
    }
    cr1=c(cr1,cr)
  }
  data3=data.frame(r=r2,cr=cr1)
  return(data3)
}

global_strength=function(re,r_max=15){
  if(!inherits(re,"mcf")){
    stop("need mcf object")
  }
  #standarided function
  nulls=attr(re,"null_value")
  mean_null=apply(nulls,1,mean,na.rm=T)
  sd_null=apply(nulls,1,sd,na.rm=T)
  adj_obs=(re$obs-mean_null)/sd_null
  sum(adj_obs[re$obs<=r_max],na.rm=T)*re$r[2]
}